Small RNA Sequencing Analysis of STZ-Injured Pancreas Reveals Novel MicroRNA and Transfer RNA-Derived RNA with Biomarker Potential for Diabetes Mellitus

MicroRNAs (miRNAs) and transfer RNA-derived small RNAs (tsRNAs) play critical roles in the regulation of different biological processes, but their underlying mechanisms in diabetes mellitus (DM) are still largely unknown. This study aimed to gain a better understanding of the functions of miRNAs and tsRNAs in the pathogenesis of DM. A high-fat diet (HFD) and streptozocin (STZ)-induced DM rat model was established. Pancreatic tissues were obtained for subsequent studies. The miRNA and tsRNA expression profiles in the DM and control groups were obtained by RNA sequencing and validated with quantitative reverse transcription-PCR (qRT-PCR). Subsequently, bioinformatics methods were used to predict target genes and the biological functions of differentially expressed miRNAs and tsRNAs. We identified 17 miRNAs and 28 tsRNAs that were significantly differentiated between the DM and control group. Subsequently, target genes were predicted for these altered miRNAs and tsRNAs, including Nalcn, Lpin2 and E2f3. These target genes were significantly enriched in localization as well as intracellular and protein binding. In addition, the results of KEGG analysis showed that the target genes were significantly enriched in the Wnt signaling pathway, insulin pathway, MAPK signaling pathway and Hippo signaling pathway. This study revealed the expression profiles of miRNAs and tsRNAs in the pancreas of a DM rat model using small RNA-Seq and predicted the target genes and associated pathways using bioinformatics analysis. Our findings provide a novel aspect in understanding the mechanisms of DM and identify potential targets for the diagnosis and treatment of DM.


Introduction
Diabetes mellitus (DM) is a chronic disease characterized by abnormal blood glucose levels and associated with short-and long-term complications [1]. According to the World Health Organization, diabetes is estimated to become the seventh leading cause of death worldwide by 2030 [2]. In order to further study the etiology, pathogenesis and effective treatment, animal models for diabetes are needed in basic experimental research. To date, some animal models have been established to study diabetes. In particular, the rat model induced by HFD and STZ is widely used. Firstly, rats are fed with a high-fat diet to induce obesity, lipid metabolism disorder and insulin resistance [3,4]. Following this, low-dose STZ is injected intravenously or intraperitoneally to damage part of the pancreatic tissue

The H&E Staining of the Pancreas in DM Model and Normal Rats
Rat pancreatic tissues were morphologically examined by H&E staining. The microscopic findings indicated that there was damage to the pancreas of DM rats. For example, the number of pancreatic cells was reduced, and the diameter of the islets became shorter. In addition, the islets were structurally disorganized, with vacuoles and swollen nuclei ( Figure 1).
In addition, the islets were structurally disorganized, with vacuoles and swollen nuclei ( Figure 1).

Profiling and Characterization of tsRNAs in Normal and Model Groups
To investigate whether the expression of tsRNAs is altered in the pancreatic tissue of diabetic rats, we analyzed the number of subtypes of tsRNA transcripts in the normal and model groups. The pie charts show the distribution of the number of each subtype of tsRNA in the normal and model groups (Figure 2A,B). In addition, Figure 2C,D shows the subtype of tsRNAs derived from different anticodons. In the Venn plot, a total of 324 expressed tsRNAs are shown in the normal and model groups, in addition to 45 and 15 specifically expressed tsRNAs in the normal and model groups, respectively ( Figure 2E). Figure 2F shows that the tsRNAs detected in our experiment were all unknown tsRNAs in the tRFdb (relational database of transferred RNA-associated fragments).

Profiling and Characterization of tsRNAs in Normal and Model Groups
To investigate whether the expression of tsRNAs is altered in the pancreatic tissue of diabetic rats, we analyzed the number of subtypes of tsRNA transcripts in the normal and model groups. The pie charts show the distribution of the number of each subtype of tsRNA in the normal and model groups (Figure 2A,B). In addition, Figure 2C,D shows the subtype of tsRNAs derived from different anticodons. In the Venn plot, a total of 324 expressed tsRNAs are shown in the normal and model groups, in addition to 45 and 15 specifically expressed tsRNAs in the normal and model groups, respectively ( Figure 2E). Figure 2F shows that the tsRNAs detected in our experiment were all unknown tsRNAs in the tRFdb (relational database of transferred RNA-associated fragments).

Differentially Expressed miRNAs and tsRNAs between Normal and Model Groups
We represented the expression changes of miRNAs and tsRNAs in pancreatic tissues of diabetic rats in the form of heat map, scatter plot and volcano plot ( Figure 3A,B). A total of 17 differentially expressed miRNAs were found in the pancreatic tissues of the model group compared with the normal group, including 9 up-and 8 down-regulated miRNAs ( Figure 3E). Among them, the most significantly up-and down-regulated miRNAs in the model group were rno-miR-19a-5p and rno-miR-3594-5p, respectively (Table 2). Furthermore, we found that a total of 11 tsRNAs were significantly up-regulated and 17 tsRNAs were significantly down-regulated in the model group compared with the normal group ( Figure 3F). Among them, tRF-Arg-ACG-007 (Fold Change ≈ 9.85) and tRF-Ser-GCT-008 (Fold Change ≈ 0.05) were the most significant up/down-regulated tsRNAs ( Table 3).

Differentially Expressed miRNAs and tsRNAs between Normal and Model Groups
We represented the expression changes of miRNAs and tsRNAs in pancreatic tissues of diabetic rats in the form of heat map, scatter plot and volcano plot ( Figure 3A,B). A total of 17 differentially expressed miRNAs were found in the pancreatic tissues of the model group compared with the normal group, including 9 up-and 8 down-regulated miRNAs ( Figure 3E). Among them, the most significantly up-and down-regulated miRNAs in the model group were rno-miR-19a-5p and rno-miR-3594-5p, respectively (Table 2). Furthermore, we found that a total of 11 tsRNAs were significantly up-regulated and 17 tsRNAs were significantly down-regulated in the model group compared with the normal group ( Figure 3F). Among them, tRF-Arg-ACG-007 (Fold Change ≈ 9.85) and tRF-Ser-GCT-008 (Fold Change ≈ 0.05) were the most significant up/down-regulated tsRNAs ( Table 3).  . The color scale is as follows: blue represents an expression level below the mean, and red represents an expression level above the mean. The colored bar top of the top panel shows the sample group, and the colored bar at the right side of the panel indicates the divisions that were performed using K-means. (C,D) The CPM values of all miRNAs and tsRNAs are plotted. The values of X and Y axes in the scatter plot are the averaged CPM values of each group (log2 scaled). miRNAs and tsRNAs above the top line (red dots, up-regulation) or below the bottom line (green dots, down-regulation) indicate more than a 1.5-fold change between the two compared groups. Gray dots indicate non-differentially expressed miRNAs and tsRNAs. (E,F) The values of X and Y axes in the volcano plot are log2 transformed fold change and −log10 transformed p-values between the two groups, respectively. Red/green circles indicate statistically significant differentially expressed miRNAs and tsRNAs with a fold change no less than 1.5 and p-value ≤ 0.05 (red: up-regulated; green: down-regulated). Gray circles indicate non-differentially expressed miRNAs and tsRNAs, with an FC and/or q-value not meeting the cut-off thresholds.

Identification of miRNA and tsRNA Target Genes
Since miRNAs and tsRNAs can function by acting on their source genes, we constructed a miRNA/tsRNA-target gene network. First, based on the prediction results, we found that 17 miRNAs and 12 tsRNAs were associated with 5186 and 2521 target genes, respectively (Supplementary Figures S1-S4). We further retained the target genes with a connectivity degree larger than three to construct the interaction network. As shown in Figure 5, a total of 35 and 19 target genes interacted with up-and down-regulated miRNAs, respectively, while a total of 5 and 17 target genes interacted with up-and down-regulated tsRNAs, respectively. These target genes, including Lpin2, Erf3 and Fkbp5, may be closely related to the role of miRNAs and tsRNAs in the development of DM.

Identification of miRNA and tsRNA Target Genes
Since miRNAs and tsRNAs can function by acting on their source genes, we constructed a miRNA/tsRNA-target gene network. First, based on the prediction results, we found that 17 miRNAs and 12 tsRNAs were associated with 5186 and 2521 target genes, respectively (Supplementary Figures S1-S4). We further retained the target genes with a connectivity degree larger than three to construct the interaction network. As shown in Figure 5, a total of 35 and 19 target genes interacted with up-and down-regulated miR-NAs, respectively, while a total of 5 and 17 target genes interacted with up-and downregulated tsRNAs, respectively. These target genes, including Lpin2, Erf3 and Fkbp5, may be closely related to the role of miRNAs and tsRNAs in the development of DM.

Enrichment Analysis for the Target Genes of Differentially Expressed miRNAs and tsRNAs
To gain additional insight into the potential mechanisms of differentially expressed miRNAs and tsRNAs, we conducted GO and KEGG pathway enrichment analyses. The results of GO analysis of differentially expressed miRNA target genes are shown in Figure  6A,B, in which the most significantly enriched Biological Process (BP) is localization (GO:

Enrichment Analysis for the Target Genes of Differentially Expressed miRNAs and tsRNAs
To gain additional insight into the potential mechanisms of differentially expressed miRNAs and tsRNAs, we conducted GO and KEGG pathway enrichment analyses. The results of GO analysis of differentially expressed miRNA target genes are shown in Figure 6A,B, in which the most significantly enriched Biological Process (BP) is localization (GO: 0051179), the Cellular Component (CC) is intracellular (GO: 0005622) and the Molecular Function (MF) is protein binding (GO: 0005515). In addition, the target genes of down-regulated expressed tsRNAs were significantly enriched in BP, CC and MF for signaling regulation (GO:0023051), intracellular (GO:0005622) and binding (GO:0005488), respectively. The target genes of up-regulated expressed tsRNAs were mainly enriched in localization (GO:0051179), intracellular (GO:0005622) and protein binding (GO:0005515). For the KEGG pathway enrichment analysis, we found that the most enriched pathways for the target genes of down-and up-regulated expressed miRNAs are leukocyte trans-endothelial migration ( Figure 7A) and choline metabolism in cancer ( Figure 7B), respectively. In addition, the most predominant enriched pathways for the target genes of down-and up-regulated expressed tsRNAs are butyrate metabolism ( Figure 7C) and proteoglycan in cancer ( Figure 7D), respectively. the KEGG pathway enrichment analysis, we found that the most enriched pathways for the target genes of down-and up-regulated expressed miRNAs are leukocyte trans-endothelial migration ( Figure 7A) and choline metabolism in cancer ( Figure 7B), respectively. In addition, the most predominant enriched pathways for the target genes of down-and up-regulated expressed tsRNAs are butyrate metabolism ( Figure 7C) and proteoglycan in cancer ( Figure 7D), respectively.

Discussion
tsRNAs are a class of Dicer-dependent noncoding RNAs that are found in various cell types [20]. Previous studies have shown that tsRNAs have a regulatory role in physiological and pathological processes [21]. Similar to the mechanism of action of miRNAs, tsRNAs can also regulate RNA stability by binding to mRNAs and thus participate in the regulation of physiopathological processes [22].

Discussion
tsRNAs are a class of Dicer-dependent noncoding RNAs that are found in various cell types [20]. Previous studies have shown that tsRNAs have a regulatory role in physiological and pathological processes [21]. Similar to the mechanism of action of miRNAs, tsRNAs can also regulate RNA stability by binding to mRNAs and thus participate in the regulation of physiopathological processes [22].
In order to investigate the expression profiles of miRNA and tsRNA and their function, we performed RNA-seq on the pancreas using a DM rat model. In this study, 17 miRNAs and 28 tsRNAs were differentially expressed between model and normal groups. Several miRNAs are involved in the pathogenesis of DM, including miR-542-5p, miR-34a-5p, miR-182, miR-217-5p, miR-448-3p and miR-384-5p [23][24][25][26][27][28][29][30]. Among these miRNAs, miR-542-5p and miR-182 could lead to potential hyperglycemia and modulate the insulin secretion by targeting FOXO1 [23,31]. miR-34a-5p may induce some expression change in PPARγ, which is a main regulator of lipid and glucose metabolism [32]. In addition, mitochondria play a crucial role in metabolic homeostasis, and alteration of mitochondrial function is a hallmark of diabetes [33]. The research by Roser Farre Garros et al. showed that upregulated miR-542-5p potentially contribute to muscle atrophy by promoting mitochondrial dysfunction [34]. miR-217-5p could affect mitochondrial function and regulate energy metabolism by targeting SIRT1 [35]. Furthermore, miR-34a-5p and miR-182 might regulate the processes of cell proliferation and apoptosis and also be involved in the modulation of beta cell survival [31,36,37]. Therefore, the roles played by miRNAs may be related to beta cell survival, lipid metabolism, glucose metabolism and mitochondrial function. In our study, miR-19a-5p and miR-3594-5p were the top up-and down-regulated miRNA in the model group as determined by RNA-Seq. However, the role of these miRNAs in the pathogenesis of DM have not been fully investigated. Further experiments will be required to address these questions.
It is known that miRNAs and tsRNAs can regulate gene expression at different levels [38,39]. To better understand the function of miRNAs and tsRNAs, we predicted the target genes of miRNAs and tsRNAs and constructed the interaction network. In the present study, a total of 5186 and 2521 mRNAs interacted with significantly altered miRNAs and tsRNAs respectively. We found that miRNAs may participate in the development of DM by regulating Nalcn, Lpin2 and Hoxa11, which were associated with insulin release, fat distribution and inflammation [40][41][42]. Among the 1719 target genes of the four up-regulated tsRNAs, five genes, including E2f3, Wdr41, Fkbp5, Ubash3a and Nr3c2, were associated with three tsRNAs. Among these genes, Ubash3a was a candidate risk factor in type 1 diabetes [43]. Wdr41 was associated with the development of T2DM [44]. The expression of E2f3 has a role in high-glucose-induced vascular endothelial injury and β-cell quiescence and proliferation [45,46]. In addition, diabetes is usually accompanied by dysregulation of the hypothalamic-pituitary-adrenal (HPA) axis, which is centrally regulated through glucocorticoid (GR) and mineralocorticoid receptors (MR). In the study of Olaf et al., Nr3c2, encoding the mineralocorticoid receptor, was up-regulated in Zucker diabetic fatty rats and implicated in the dysregulation of the HPA axis [47]. In a study of 20 T2DM subjects and 20 non-diabetic subjects, Cherno et al. suggested that Fkbp51 might be a key factor in glucocorticoid-induced insulin resistance [48]. On the other hand, among the 1351 target genes of the eight down-regulated tsRNAs, Irs2 was associated with three tsRNAs. Irs2 is a major component of the insulin/insulin-like growth factor-1 signaling pathway. The disruption of Irs2 impairs both peripheral insulin signaling and pancreatic β-cell function and leads to life-threatening T2DM [49]. In summary, the differentially expressed miRNAs and tsRNAs may participate in the development of DM through regulating the above key genes.
To further reveal the metabolic pathways associated with the target genes of differentially expressed miRNAs and tsRNAs, we performed KEGG pathway enrichment analysis. The results showed that the target genes of miRNAs are enriched in the Wnt signaling pathway, and the target genes of tsRNAs are enriched in the insulin signaling pathway. Wnt proteins, a protein family composed of secreted glycoproteins, are able to regulate a variety of developmental processes [50]. Previous research found that Wnt signaling regulates pancreatic islet proliferation and β-cell proliferation and is crucial for endocrine development [51,52]. In addition, the MAPK signaling pathway and Hippo signaling pathway were associated with both miRNAs and tsRNAs. The MAPK signaling pathway can mediate expression of vascular endothelial growth factor (VEGF), which is one of the major regulatory molecules in diabetes [53]. The MAPK pathway also plays an important role in the inflammation process [54]. Under diabetic conditions, the MAPK signaling pathway is activated and causes the expression of related inflammatory factors [55]. It is closely related to Ras. Hippo signaling is an evolutionarily conserved pathway that critically regulates development and homeostasis of various tissues. This pathway has been proven to play important roles in pancreas development and the regulation of β-cell survival, proliferation and regeneration [56]. Therefore, based on the results of KEGG analysis, we speculated that the altered miRNAs and tsRNAs may influence the development of DM through functioning in pathways associated with inflammation, glycan metabolism and endocrine development.
There are a few limitations of the present study. Firstly, only one timepoint after modeling was selected for study. Therefore, we could not observe the alteration of miRNAs and tsRNAs in the development of DM. Secondly, we selected the pancreas as the study subject. However, due to the complexity of DM pathogenesis, other tissues should be evaluated in future studies. Last but not least, we did not investigate the definitive molecular mechanisms of the differentially expressed miRNAs and tsRNAs, which will be examined in our future projects.
Taken together, this study revealed the expression profiles of miRNAs and tsRNAs in the pancreas using a DM rat model and small RNA-Seq and predicted the target genes and the pathways involved using bioinformatics analyses. Our findings suggest the underlying mechanisms of DM as well as potential targets to the diagnosis and treatment of DM.

Animal Experiment and Tissue Extraction
All experimental protocols were approved by the Animal Care and Management Committee of Beijing University of Chinese Medicine. Male Wistar rats (8 weeks old, Beijing Life River Laboratory Animal Technology Co., Ltd., Beijing, China) were used in this study. All rats were housed in plastic cages at 22 ± 2 • C with a 12-h light/dark cycle. After one week of acclimatization feeding, 20 rats were randomly divided into two groups: the normal group and model group. The rats of the normal group were fed a standard diet (AIN-96G diet, Sibeifu Bioscience Co., Ltd., Beijing, China), and the rats of the model group were fed a high-fat diet (HFD, D12451, 45% fat Kcal%, Sibeifu Bioscience Co., Ltd., Beijing, China) for six weeks. After 12 h of fasting, HFD-fed rats were injected with STZ 30 mg/kg to induce the DM model. Control rats were then injected with an equal volume of sterile sodium citrate buffer. After one week, HFD-fed rats were considered to have DM if fasting blood glucose (FBG) was ≥11.1 mmol/L. Finally, all rats were anesthetized by intraperitoneal injection of pentobarbital sodium anesthesia at 40 mg/kg. After that, rats were sacrificed by exsanguination of the femoral artery. Blood samples and pancreatic tissues were obtained for further experiments.

Histological Examination of Pancreatic Tissues
Pancreatic tissues were fixed with formaldehyde and embedded in paraffin wax. After sectioning, the slices were stained with hematoxylin and eosin (H&E) according to established procedures. Finally, histopathological changes in the pancreatic tissue were observed with an optical microscope (Olympus, Tokyo, Japan).

RNA Isolation and Small RNA-Seq Anaysis
RNA was isolated from the pancreatic tissues of diabetic and normal rats using the Trizol RNA extraction kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Each RNA sample was subsequently checked for integrity and concentration using agarose gel electrophoresis and the Nanodrop TM instrument (Thermo Fisher Scientific, Waltham, MA, USA). Total RNA from each sample was ligated to the 3 and 5 small RNA aptamers. The cDNA was then synthesized and amplified using Illumina proprietary RT primers and amplification primers. Afterwards, 134-160 bp PCR amplicons were extracted and purified from PAGE gels. Finally, total RNA was quantified using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing was performed on the Illumina NextSeq 500 system using the NextSeq 500/550 V2 kit (#FC-404-2005, Illumina, San Diego, CA, USA) according to the manufacturer's instructions.

Data Collection and Analysis
Raw sequencing data generated in the form of Illumina NextSeq 500 was trimmed with cutadapt for sequencing reads with a 5-and 3-adaptor and discarded reads (length <14 nt or length >40 nt) after Illumina quality control. The reads were then compared using bowtie software and miRDeep, respectively, to allow only 1 mismatch to the precursor tRNA sequence and the precursor tRNA sequence. Based on the statistical analysis of the comparison (mapping rate, read length, fragment sequence deviation), we determined whether the results could be used for subsequent data analysis. Finally, the expression profiles of tsRNAs and miRNAs in the normal and model group were obtained.

Validation of Small RNA-Seq Analysis Using qRT-PCR
The sequencing results of selected tsRNAs and miRNAs were validated by qRT-PCR. RNA samples for tsRNAs were firstly pretreated using the rtStar TM tRF&tiRNA Pretreatment Kit (Arraystar, Rockville, MD, USA). Then, cDNA synthesis was performed with the rtStar TM First Strand cDNA Synthesis Kit (Arraystar). Afterwards, the samples were detected using the following thermal cycling parameters on the ViiA 7 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA): initial activation for 10 min at 95 • C, followed by 40 cycles of denaturation for 10 s at 95 • C, and annealing and extension for 60 s at 60 • C. RNA samples of miRNAs were reverse-transcribed with MMLV reverse transcriptase (Epicentre Technologies, Madison, WI, USA) and amplified in the Gene Amp PCR System 9700 (Applied Biosystems). Finally, the expression levels of tsRNA and miRNA were calculated by the 2 −∆∆Ct method, and U6 was used as the standard. A list of primers is given in Table 4. 4.6. Differential tsRNA and miRNA Target Gene Prediction and Functional Prediction miRanda and Targetscan were utilized to further explore the target genes of the differentially expressed tsRNAs and miRNAs. Gene Ontology (GO) (http://www.geneontology.org/, accessed on 17 September 2019) analysis aimed to obtain the relationship between the target genes and functional terms. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/, accessed on 17 September 2019) was used for pathway analysis, to find the associated pathways affected by the target genes. For both GO and pathway analysis, p < 0.05 was considered as statistically significant.

Statistical Analysis
The read counts in this study were normalized to counts per million mapped reads (CPM). The 'EdgeR' R package was used for differential expression analysis. Differentially expressed miRNAs and tsRNAs were identified through fold-change screening (fold change ≥1.5 and p < 0.05). SPSS version 22.0 (SPSS Inc., Chicago, IL, USA) was used for the statistics. Data were expressed as mean ± SEM and compared using the unpaired Student's t-test. p < 0.05 was considered statistically significant.